S1 Text: Modelling results
Comparative analysis of dengue and Zika outbreaks reveals differences by setting and virus
Sebastian Funk, Adam J. Kucharski, Anton Camacho, Rosalind M. Eggo, Laith Yakob, Lawrence Murray, W. John Edmunds
Introduction
This document contains the R code necessary to reproduce the results in “Comparative analysis of dengue and Zika outbreaks reveals differences by setting and virus.” All code and data used here are in a github repository. The data sets for the three outbreaks are contained in the data subdirectory.
Required packages
cran_packages <- c("cowplot", "dplyr", "tidyr", "stringi", "msm", "rbi")
github_packages <- c("sbfnk/RBi.helpers")for (package in cran_packages) {
install.packages(package)
}
library("devtools")
for (package in github_packages) {
install_github(package)
}Generate MCMC chains
Generating MCMC chains requires a working installation of libbi v1.2.0. The model code is in vbd.bi, which is run by the R script dengue_zika_mcmc.r.
Edit the lines starting code_dir <- and output_dir <- in the script ‘dengue_zika_mcmc.r’ from the R/ directory in the github repository to point to the correct directories, and run it using
Rscript dengue_zika_mcmc.r --nsamples 10000000 --pre_samples 10000 --sample-prior --sample-observations --thin 1000 --fix-natural-history
Rscript dengue_zika_mcmc.r --nsamples 10000000 --pre_samples 10000 --sample-prior --sample-observations --thin 1000 --fix-natural-history --shorter-mosquito-lifespanThe first line runs with \(D_\mathrm{life,M}=2\,\mathrm{weeks}\), and the second one with \(D_\mathrm{life,M}=1\,\mathrm{week}\)
Analyse MCMC chains
Set the code_dir and output_dir to the value used in sample.posterior.r
output_dir <- "insert your directory here"
code_dir <- "insert your directory here"Load required packages
for (package in c(cran_packages, github_packages)) {
library(tolower(sub("^.*/", "", package)), character.only = TRUE)
}Load script from github repository (add path if necessary)
source(paste0(code_dir, "/R/", "analyse_traces.r"))Get MCMC chains of the two libbi runs.
libbi_results <- analyse_traces(c("vbd_fnh", "vbd_fnh_shorter"), output_dir)Calculate DIC for models with different mosquito life times (1 week vs. 2 weeks vs. both with equal probability)
dic <- sapply(libbi_results, function(x) {
compute_DIC(x[["trace"]][["posterior"]])
})
dic## vbd_fnh vbd_fnh_shorter vbd_fnh_all
## 366.7032 365.8579 366.2954Plots
Prepare data
ts <- list()
analyses <- data.frame(setting = c("yap", "yap", "fais"), disease = c("dengue",
"zika", "dengue"))
for (i in 1:nrow(analyses)) {
this_setting <- analyses[i, "setting"]
this_disease <- analyses[i, "disease"]
this_filename <- paste(code_dir, "data", paste(this_setting, this_disease,
"data.rds", sep = "_"), sep = "/")
this_ts <- readRDS(this_filename) %>% mutate(setting = this_setting, disease = this_disease,
week = floor(nr/7))
ts <- c(ts, list(this_ts))
}
ordered_obs_id_levels <- c("yap_zika", "yap_dengue", "fais_dengue")
data_labels <- ordered_obs_id_levels
data_labels <- sub("^(.*)_(.*)$", "\\2 \\1", data_labels)
data_labels <- sub(" ", " in ", stri_trans_totitle(data_labels))
names(data_labels) <- ordered_obs_id_levels
data <- bind_rows(ts) %>% group_by(week, setting, disease) %>% summarize(value = sum(value),
onset_date = min(onset_date)) %>% ungroup() %>% mutate(obs_id = factor(paste(setting,
disease, sep = "_"), levels = ordered_obs_id_levels, labels = data_labels)) %>%
arrange(week, obs_id) %>% select(week, obs_id, value, onset_date) %>% rename(time = week) %>%
mutate(state = "Cases")
first_obs <- data %>% group_by(obs_id) %>% filter(value > 0) %>% slice(which.min(time)) %>%
select(time, obs_id) %>% rename(first_obs = time)
last_obs <- data %>% group_by(obs_id) %>% filter(value > 0) %>% slice(which.max(time)) %>%
select(time, obs_id) %>% rename(last_obs = time)
data <- data %>% left_join(first_obs, by = "obs_id") %>% left_join(last_obs,
by = "obs_id") %>% filter(time >= first_obs & time <= last_obs)Fits
libbi_results[["vbd_fnh_all"]][["trace"]][["posterior"]][["Cases"]] <- libbi_results[["vbd_fnh_all"]][["trace"]][["posterior"]][["Cases"]] %>%
mutate(obs_id = factor(obs_id, levels = ordered_obs_id_levels, labels = data_labels))
temp_plot <- plot_libbi(read = libbi_results[["vbd_fnh_all"]][["trace"]][["posterior"]]["Cases"],
model = libbi_results[["vbd_fnh_all"]][["model"]], data = data %>% filter(value >
0), density_args = list(adjust = 2), extra.aes = list(group = "obs_id"),
data.colour = "black", states = "Cases", trend = "mean", plot = FALSE, limit.to.data = TRUE,
quantiles = c(0.5, 0.72, 0.95))
obs_states <- temp_plot$data$states %>% inner_join(data %>% select(time, obs_id,
onset_date), by = c("time", "obs_id"))
p_obs <- ggplot(obs_states, aes(x = onset_date)) + geom_point(data = data, mapping = aes(y = value)) +
facet_wrap(~obs_id, scales = "free") + scale_x_date("Week", labels = scales::date_format("%e %b %Y")) +
theme(axis.text.x = element_text(angle = 45, hjust = 1)) + facet_wrap(~obs_id,
scales = "free") + scale_y_continuous("Disease incidence") + geom_line(aes(y = value)) +
geom_ribbon(aes(ymin = min.1, ymax = max.1), alpha = 0.5) + geom_ribbon(aes(ymin = min.2,
ymax = max.2), alpha = 0.25) + geom_ribbon(aes(ymin = min.3, ymax = max.3),
alpha = 0.125)
p_obs
Other plots
Define plot labels
labels <- c(p_d_inc_h = "italic(D)[plain(inc,H)]", p_d_inc_m = "italic(D)[plain(inc,M)]",
p_d_inf_h = "italic(D)[plain(inf,H)]", p_lm = "log[10](italic(m))", p_initial_susceptible_yap = "italic(q)",
p_rep = "italic(r)", p_b_h = "italic(b)[H]", p_b_m = "italic(b)[M]", p_t_start = "italic(t[0])",
p_pop_yap = "italic(rho)", p_red_foi_yap = "italic(sigma)", p_p_patch_yap = "italic(phi)",
R0 = "italic(R)[H %->% H]", GI = "italic(G)", zika = "Zika", yap = "Yap",
fais = "Fais")Marginal posterior densities
p <- plot_libbi(read = libbi_results[["vbd_fnh_all"]][["trace"]][["posterior"]],
prior = libbi_results[["vbd_fnh_all"]][["trace"]][["prior"]], model = libbi_results[["vbd_fnh_all"]][["model"]],
density_args = list(bins = 20, alpha = 0.5, color = "black"), densities = "histogram",
extra.aes = list(color = "disease", linetype = "setting"), trend = "median",
plot = FALSE, quantiles = c(0.5, 0.95), labels = labels, brewer.palette = "Set1")
p$densities
Human-to-human reproduction number vs generation interval
temp_plot <- plot_libbi(read = libbi_results[["vbd_fnh_all"]][["trace"]][["posterior"]],
prior = libbi_results[["vbd_fnh_all"]][["trace"]][["prior"]], model = libbi_results[["vbd_fnh_all"]][["model"]],
extra.aes = list(color = "disease", linetype = "setting", group = "p_d_life_m"),
plot = FALSE, labels = labels, states = c(), params = c("R0", "GI"), noises = c())
r0gi <- temp_plot$data$params %>% filter(distribution == "posterior") %>% mutate(obs_id = tolower(paste(setting,
disease, sep = "_"))) %>% filter(obs_id != "fais_zika") %>% mutate(obs_id = factor(obs_id,
levels = ordered_obs_id_levels, labels = data_labels)) %>% spread(parameter,
value)
cross_sections <- data.frame(GI = c(3, 4))
colnames(cross_sections) <- "italic(G)"
p_r0vgi <- ggplot(r0gi, aes(x = `italic(G)`)) + geom_jitter(aes(y = `italic(R)[H %->% H]`,
color = factor(p_d_life_m))) + facet_grid(~obs_id) + scale_x_continuous("Equilibrium generation interval (weeks)") +
scale_y_continuous(expression(italic(R)[H %->% H])) + scale_color_brewer("Mosquito life span",
palette = "Dark2", labels = c("1 week", "2 weeks")) + theme(legend.position = "top") +
facet_wrap(~obs_id, scales = "free") + geom_rect(data = cross_sections,
aes(xmin = `italic(G)` - 0.1, xmax = `italic(G)` + 0.1), ymin = -Inf, ymax = Inf,
alpha = 0.2, lwd = 3)
p_r0vgi
Traces
p$traces
Other models
These are the alternative models considered and mentioned in the manuscript.
Reduced population size in Yap
Generate MCMC chains
In this model, the population in Yap is reduced by a factor \(rho\).
Fit the model (remember to update code_dir and data_dir in dengue_zika_mcmc.r):
Rscript dengue_zika_mcmc.r --nsamples 10000000 --pre_samples 10000 --sample-prior --sample-observations --thin 1000 --fix-natural-history --sero --pop
Rscript dengue_zika_mcmc.r --nsamples 10000000 --pre_samples 10000 --sample-prior --sample-observations --thin 1000 --fix-natural-history --shorter-mosquito-lifespan --sero --popGet MCMC chains of the two libbi runs.
libbi_results <- analyse_traces(c("vbd_sero_pop_fnh", "vbd_sero_pop_fnh_shorter"),
output_dir)Fits
libbi_results[["vbd_sero_pop_fnh_all"]][["trace"]][["posterior"]][["Cases"]] <- libbi_results[["vbd_sero_pop_fnh_all"]][["trace"]][["posterior"]][["Cases"]] %>%
mutate(obs_id = factor(obs_id, levels = ordered_obs_id_levels, labels = data_labels))
temp_plot <- plot_libbi(read = libbi_results[["vbd_sero_pop_fnh_all"]][["trace"]][["posterior"]]["Cases"],
model = libbi_results[["vbd_sero_pop_fnh_all"]][["model"]], data = data %>%
filter(value > 0), density_args = list(adjust = 2), extra.aes = list(group = "obs_id"),
data.colour = "black", states = "Cases", trend = "mean", plot = FALSE, limit.to.data = TRUE,
quantiles = c(0.5, 0.72, 0.95))
obs_states <- temp_plot$data$states %>% inner_join(data %>% select(time, obs_id,
onset_date), by = c("time", "obs_id"))
p_obs <- ggplot(obs_states, aes(x = onset_date)) + geom_point(data = data, mapping = aes(y = value)) +
facet_wrap(~obs_id, scales = "free") + scale_x_date("Week", labels = scales::date_format("%e %b %Y")) +
theme(axis.text.x = element_text(angle = 45, hjust = 1)) + facet_wrap(~obs_id,
scales = "free") + scale_y_continuous("Disease incidence") + geom_line(aes(y = value)) +
geom_ribbon(aes(ymin = min.1, ymax = max.1), alpha = 0.5) + geom_ribbon(aes(ymin = min.2,
ymax = max.2), alpha = 0.25) + geom_ribbon(aes(ymin = min.3, ymax = max.3),
alpha = 0.125)
p_obs
Marginal posterior densities
p <- plot_libbi(read = libbi_results[["vbd_sero_pop_fnh_all"]][["trace"]][["posterior"]],
prior = libbi_results[["vbd_sero_pop_fnh_all"]][["trace"]][["prior"]], model = libbi_results[["vbd_sero_pop_fnh_all"]][["model"]],
density_args = list(bins = 20, alpha = 0.5, color = "black"), densities = "histogram",
extra.aes = list(color = "disease", linetype = "setting"), trend = "median",
plot = FALSE, quantiles = c(0.5, 0.95), labels = labels, brewer.palette = "Set1")
p$densities
Estimated final size for Zika in Yap (proportion infected in the outbreak)
n_infected <- libbi_results[["vbd_sero_pop_fnh_all"]][["trace"]][["posterior"]][["final_size"]] %>%
filter(obs_id == "yap_zika") %>% .$value %>% quantile(c(0.025, 0.5, 0.975))
N <- libbi_results[["vbd_sero_pop_fnh_all"]][["trace"]][["posterior"]][["p_N_h"]] %>%
filter(setting == "yap") %>% .$value %>% mean
n_infected/N## 2.5% 50% 97.5%
## 0.6820483 0.7246585 0.7699886Reduced susceptibility against Zika in Yap
Generate MCMC chains
In this model, the population susceptible to Zika in Yap is reduced by a factor \(r\).
Fit the model (remember to update code_dir and data_dir in dengue_zika_mcmc.r):
Rscript dengue_zika_mcmc.r --nsamples 10000000 --pre_samples 10000 --sample-prior --sample-observations --thin 1000 --fix-natural-history --sero
Rscript dengue_zika_mcmc.r --nsamples 10000000 --pre_samples 10000 --sample-prior --sample-observations --thin 1000 --fix-natural-history --shorter-mosquito-lifespan --seroGet MCMC chains of the two libbi runs
libbi_results <- analyse_traces(c("vbd_sero_fnh", "vbd_sero_fnh_shorter"), output_dir)Fits
libbi_results[["vbd_sero_fnh_all"]][["trace"]][["posterior"]][["Cases"]] <- libbi_results[["vbd_sero_fnh_all"]][["trace"]][["posterior"]][["Cases"]] %>%
mutate(obs_id = factor(obs_id, levels = ordered_obs_id_levels, labels = data_labels))
temp_plot <- plot_libbi(read = libbi_results[["vbd_sero_fnh_all"]][["trace"]][["posterior"]]["Cases"],
model = libbi_results[["vbd_sero_fnh_all"]][["model"]], data = data %>%
filter(value > 0), density_args = list(adjust = 2), extra.aes = list(group = "obs_id"),
data.colour = "black", states = "Cases", trend = "mean", plot = FALSE, limit.to.data = TRUE,
quantiles = c(0.5, 0.72, 0.95))
obs_states <- temp_plot$data$states %>% inner_join(data %>% select(time, obs_id,
onset_date), by = c("time", "obs_id"))
p_obs <- ggplot(obs_states, aes(x = onset_date)) + geom_point(data = data, mapping = aes(y = value)) +
facet_wrap(~obs_id, scales = "free") + scale_x_date("Week", labels = scales::date_format("%e %b %Y")) +
theme(axis.text.x = element_text(angle = 45, hjust = 1)) + facet_wrap(~obs_id,
scales = "free") + scale_y_continuous("Disease incidence") + geom_line(aes(y = value)) +
geom_ribbon(aes(ymin = min.1, ymax = max.1), alpha = 0.5) + geom_ribbon(aes(ymin = min.2,
ymax = max.2), alpha = 0.25) + geom_ribbon(aes(ymin = min.3, ymax = max.3),
alpha = 0.125)
p_obs
Marginal posterior densities
p <- plot_libbi(read = libbi_results[["vbd_sero_fnh_all"]][["trace"]][["posterior"]],
prior = libbi_results[["vbd_sero_fnh_all"]][["trace"]][["prior"]], model = libbi_results[["vbd_sero_fnh_all"]][["model"]],
density_args = list(bins = 20, alpha = 0.5, color = "black"), densities = "histogram",
extra.aes = list(color = "disease", linetype = "setting"), trend = "median",
plot = FALSE, quantiles = c(0.5, 0.95), labels = labels, brewer.palette = "Set1")
p$densities
Estimated final size for Zika in Yap (proportion infected in the outbreak)
n_infected <- libbi_results[["vbd_sero_fnh_all"]][["trace"]][["posterior"]][["final_size"]] %>%
filter(obs_id == "yap_zika") %>% .$value %>% quantile(c(0.025, 0.5, 0.975))
N <- libbi_results[["vbd_sero_fnh_all"]][["trace"]][["posterior"]][["p_N_h"]] %>%
filter(setting == "yap") %>% .$value %>% mean
n_infected/N## 2.5% 50% 97.5%
## 0.6844271 0.7296974 0.7736204Two-patch model
Fit the model (remember to update code_dir and data_dir in dengue_zika_mcmc.r):
Rscript dengue_zika_mcmc.r --nsamples 10000000 --pre_samples 10000 --sample-prior --sample-observations --thin 1000 --fix-natural-history --patch
Rscript dengue_zika_mcmc.r --nsamples 10000000 --pre_samples 10000 --sample-prior --sample-observations --thin 1000 --fix-natural-history --shorter-mosquito-lifespan --patchGet MCMC chains of the two libbi runs.
libbi_results <- analyse_traces(c("vbd_patch_fnh", "vbd_patch_fnh_shorter"),
output_dir)Fits
libbi_results[["vbd_patch_fnh_all"]][["trace"]][["posterior"]][["Cases"]] <- libbi_results[["vbd_patch_fnh_all"]][["trace"]][["posterior"]][["Cases"]] %>%
mutate(obs_id = factor(obs_id, levels = ordered_obs_id_levels, labels = data_labels))
temp_plot <- plot_libbi(read = libbi_results[["vbd_patch_fnh_all"]][["trace"]][["posterior"]]["Cases"],
model = libbi_results[["vbd_patch_fnh_all"]][["model"]], data = data %>%
filter(value > 0), density_args = list(adjust = 2), extra.aes = list(group = "obs_id"),
data.colour = "black", states = "Cases", trend = "mean", plot = FALSE, limit.to.data = TRUE,
quantiles = c(0.5, 0.72, 0.95))
obs_states <- temp_plot$data$states %>% inner_join(data %>% select(time, obs_id,
onset_date), by = c("time", "obs_id"))
p_obs <- ggplot(obs_states, aes(x = onset_date)) + geom_point(data = data, mapping = aes(y = value)) +
facet_wrap(~obs_id, scales = "free") + scale_x_date("Week", labels = scales::date_format("%e %b %Y")) +
theme(axis.text.x = element_text(angle = 45, hjust = 1)) + facet_wrap(~obs_id,
scales = "free") + scale_y_continuous("Disease incidence") + geom_line(aes(y = value)) +
geom_ribbon(aes(ymin = min.1, ymax = max.1), alpha = 0.5) + geom_ribbon(aes(ymin = min.2,
ymax = max.2), alpha = 0.25) + geom_ribbon(aes(ymin = min.3, ymax = max.3),
alpha = 0.125)
p_obs
Marginal posterior densities
p <- plot_libbi(read = libbi_results[["vbd_patch_fnh_all"]][["trace"]][["posterior"]],
prior = libbi_results[["vbd_patch_fnh_all"]][["trace"]][["prior"]], model = libbi_results[["vbd_patch_fnh_all"]][["model"]],
density_args = list(bins = 20, alpha = 0.5, color = "black"), densities = "histogram",
extra.aes = list(color = "disease", linetype = "setting"), trend = "median",
plot = FALSE, quantiles = c(0.5, 0.95), labels = labels, brewer.palette = "Set1")
p$densities
Fitting the Zika outbreak in Yap in isolation
Fit the model (remember to update code_dir and data_dir in dengue_zika_mcmc.r):
Rscript dengue_zika_mcmc.r --nsamples 10000000 --pre_samples 10000 --sample-prior --sample-observations --thin 1000 --fix-natural-history --patch --disease zika --setting yap
Rscript dengue_zika_mcmc.r --nsamples 10000000 --pre_samples 10000 --sample-prior --sample-observations --thin 1000 --fix-natural-history --shorter-mosquito-lifespan --patch --disease zika --setting yapGet MCMC chains of the two libbi runs.
libbi_results <- analyse_traces(c("vbd_patch_fnh_yap_zika", "vbd_patch_fnh_shorter_yap_zika"),
output_dir)Fit
libbi_results[["vbd_patch_fnh_yap_zika_all"]][["trace"]][["posterior"]][["Cases"]] <- libbi_results[["vbd_patch_fnh_yap_zika_all"]][["trace"]][["posterior"]][["Cases"]] %>%
mutate(obs_id = factor(obs_id, levels = ordered_obs_id_levels, labels = data_labels))
temp_plot <- plot_libbi(read = libbi_results[["vbd_patch_fnh_yap_zika_all"]][["trace"]][["posterior"]]["Cases"],
model = libbi_results[["vbd_patch_fnh_yap_zika_all"]][["model"]], data = data,
density_args = list(adjust = 2), extra.aes = list(group = "obs_id"), data.colour = "black",
states = "Cases", trend = "mean", plot = FALSE, limit.to.data = TRUE, quantiles = c(0.5,
0.72, 0.95))
obs_states <- temp_plot$data$states %>% inner_join(data %>% select(time, obs_id,
onset_date), by = c("time", "obs_id"))
p_obs <- ggplot(obs_states, aes(x = onset_date)) + geom_point(data = data %>%
filter(obs_id == "Zika in Yap"), mapping = aes(y = value)) + facet_wrap(~obs_id,
scales = "free") + scale_x_date("Week", labels = scales::date_format("%e %b %Y")) +
theme(axis.text.x = element_text(angle = 45, hjust = 1)) + facet_wrap(~obs_id,
scales = "free") + scale_y_continuous("Disease incidence") + geom_line(aes(y = value)) +
geom_ribbon(aes(ymin = min.1, ymax = max.1), alpha = 0.5) + geom_ribbon(aes(ymin = min.2,
ymax = max.2), alpha = 0.25) + geom_ribbon(aes(ymin = min.3, ymax = max.3),
alpha = 0.125)
p_obs
Marginal posterior densities
p <- plot_libbi(read = libbi_results[["vbd_patch_fnh_yap_zika_all"]][["trace"]][["posterior"]],
prior = libbi_results[["vbd_patch_fnh_yap_zika_all"]][["trace"]][["prior"]],
model = libbi_results[["vbd_patch_fnh_yap_zika_all"]][["model"]], density_args = list(bins = 20,
alpha = 0.5, color = "black"), densities = "histogram", extra.aes = list(color = "disease",
linetype = "setting"), trend = "median", plot = FALSE, quantiles = c(0.5,
0.95), labels = labels, brewer.palette = "Set1")
p$densities
Human-to-human reproduction number vs generation interval
temp_plot <- plot_libbi(read = libbi_results[["vbd_patch_fnh_yap_zika_all"]][["trace"]][["posterior"]],
prior = libbi_results[["vbd_patch_fnh_yap_zika_all"]][["trace"]][["prior"]],
model = libbi_results[["vbd_patch_fnh_yap_zika_all"]][["model"]], extra.aes = list(color = "disease",
linetype = "setting", group = "p_d_life_m"), plot = FALSE, labels = labels,
states = c(), params = c("R0", "GI"), noises = c())
r0gi <- temp_plot$data$params %>% filter(distribution == "posterior") %>% mutate(obs_id = tolower(paste(setting,
disease, sep = "_"))) %>% filter(obs_id != "fais_zika") %>% mutate(obs_id = factor(obs_id,
levels = ordered_obs_id_levels, labels = data_labels)) %>% spread(parameter,
value)
p_r0vgi <- ggplot(r0gi, aes(x = `italic(G)`)) + geom_jitter(aes(y = `italic(R)[H %->% H]`,
color = factor(p_d_life_m))) + facet_grid(~obs_id) + scale_x_continuous("Equilibrium generation interval (weeks)") +
scale_y_continuous(expression(italic(R)[H %->% H])) + scale_color_brewer("Mosquito life span",
palette = "Dark2", labels = c("1 week", "2 weeks")) + theme(legend.position = "top") +
facet_wrap(~obs_id, scales = "free") + geom_rect(data = cross_sections,
aes(xmin = `italic(G)` - 0.1, xmax = `italic(G)` + 0.1), ymin = -Inf, ymax = Inf,
alpha = 0.2, lwd = 3)
p_r0vgi